function success = GetVolume(filePath, handles, i,varargin)
% getVolume retrieve volume information alongside the associated transform
% matrix, and converts the matrix to U8 to preserve memory.
% the first optional argument is the talairach data
plotPanel=handles.plotPanel;
% Read volume and header information
if exist(filePath,'file')
    [pathstr, name, ext] = fileparts(filePath);
        if strcmp(ext,'.img')
            fixASAHeader(fullfile(pathstr,[name '.hdr']));
        end
        header=spm_vol(filePath);
        V.vol=spm_read_vols(header);
        V.mat = header.mat;
else
    error(['Unfound file.']);
end
%V.vol = arrayfun(@ZeroToNaN,V.vol)
function y=ZeroToNaN(x)
    if x==0
        y=NaN;
    else
        y=x;
    end
end

% If coregistration matrix was computed, apply it to the transform matrix.
if exist(fullfile(pathstr,'coreg.dat'))
    coreg = load(fullfile(pathstr,'coreg.dat'));
    V.mat = coreg * V.mat;
end

% Check if not too far from the centre of the volume. If so, put it in the middle of
% the volume (neceassary for the coregistration).
T = header.mat(1:3,4);
tempMat = header.mat;
tempMat(1:3,4) = zeros(1,3);
center = -(tempMat*[(header.dim/2)';-1]);
center = center(1:3);
diff = abs(center-T./center);
if min(diff)<0.5
    header.mat(1:3,4) = zeros(1,3);
    header.mat(:,4) = -(header.mat*[(header.dim/2)';-1]);
    spm_write_vol(header,V.vol);
    V.mat = header.mat;
end
theta = get(handles.rotation,'userdata');
if isempty(theta)
    theta = 0;
end
V_rotated = makehgtform('xrotate',theta) * V.mat;
% Rescale to 0-255 // 0 reserved for strict 0, later displayed as
% transparent.
if i~=8 && i~=9
    minValue = min(min(min(V.vol)));
    maxValue = max(max(max(V.vol)));
    V.minmax = [minValue maxValue];
    N = size(V.vol);
    N = N(1);
    for k = 1:N
        V.vol(k,:,:) = (V.vol(k,:,:)~=0).*(255*(V.vol(k,:,:)-minValue)/(maxValue-minValue)+0);
    end
    V.vol = uint8(V.vol);
else
    V.vol = uint16(V.vol);
end
% if i==4
%     threshold = 0.2;
%      for k = 1:N
%         V.vol(k,:,:) = uint8(( (V.vol(k,:,:)<255*threshold) + (V.vol(k,:,:)>255*(1-threshold)))).*(V.vol(k,:,:));
%     end
% end
V.colormap = gray(256);
if i== 8 || i==9
    V.colormap=prism(2000);
elseif i==6
    V.colormap = jet(256);
end
V.alpha = 1;
if i==7
    V.alpha=0;
end
V.brightness = 0;
V.contrast = 0;
V.displayVolume = uint8(zeros(2,2,2));

% Store infos in relevant stack structure fields
stack = get(plotPanel,'userdata');
if length(stack) < i
    if i~=1
        for j=1:i-1 % Fill stack with empty
            if length(stack)<j
                stack = [stack struct('vol',uint8(zeros(size(V.vol))),'mat',makehgtform(),...
                                      'colormap',gray(256),'alpha',0,'brightness',0,'contrast',0,...
                                      'displayVolume',uint8(zeros(2,2,2)),'minmax',[0 0])];
            end
        end        
    end
    stack = [stack V];
else
    stack(i).vol = V.vol;
    stack(i).mat = V.mat;
    stack(i).colormap = V.colormap;
    stack(i).minmax = V.minmax;
end

 
if cuda  
     corners = zeros(4,8);
     volSize = size(V.vol);
        for k1 = 0:1
            for k2 = 0:1
                for k3 = 0:1
                    corners(:,k1*4+k2*2+k3+1) = V.mat * [k1 * volSize(1) ; k2 * volSize(2) ; k3 * volSize(3) ; 1];
                end
            end
        end
        MIN = min(corners');
        MAX = max(corners');
        MAX-MIN;
    dispSize = int32(floor((MAX - MIN) * 2));
   if i~=8 && i~=9
    calllib('cudaicem','ReadVolume',i,V.alpha,V.brightness,V.contrast,single(V.colormap),...
            single(V_rotated) ,size(V.vol),V.vol);
   else
       calllib('cudaicem','ReadTalairachVolume',i==9,V.alpha,V.brightness,V.contrast,single(V.colormap),...
            single(V_rotated) ,size(V.vol),V.vol,varargin{1}');
   end
end
set(plotPanel,'userdata',stack);
success = true;
end
    
